*! version 1.1.0  10/21/93   STB-17 ssi5.2
program define ridder 
	version 3.1
	parse "`*'", parse(" X=")

	local toload "exp1"
	while "`1'"!="=" & "`1'"!="returns" & "`1'"!="" {
		if "`1'"=="X" {
			if "`toload'"=="exp2" {
				di in red "X may be specifed only once"
				exit 198
			}
			local toload "exp2"
		}
		else {
			local `toload' "``toload'' `1'"
		}
		mac shift
	}
	if "`toload'"!="exp2" {
		di in red "X must be specified"
		exit 198
	}
	if "`1'"=="" { error 198 }
	if "`1'"=="returns" {
		if "`2'"=="macro" {
			local type "`3'"
			local rtnval
		}
		else if "`2'"=="exp" {
			local rtnval "`3'"
			local type
		}
		else	error 198
		mac shift 3
		local progttl "`exp1' X `exp2'"
		local inline "*"
	}
	else {
		local progttl = substr("`exp2'",2,.)
		local progttl "`exp1'X`progttl'"
		local callout "*"
	}
	if "`1'" != "=" { 
		di in red "=exp required"
		exit 198
	}
	mac shift
	local toload "mval"
	while "`1'"!="" {
		if "`1'"=="to" {
			local toload "mhi"
		}
		else if "`1'"=="from" {
			local toload "mlo"
		}
		else if "`1'"=="tol" {
			local toload "mtol"
		}
		else if "`1'"=="ltol" {
			local toload "mltol"
		}
		else {
			local `toload' "``toload''`1'"
		}
		mac shift
	}
	if "`mval'"=="" { 
		di in red "= must be followed by expression"
		exit 198
	}
	if "`mhi'"=="" { 
		di in red "to specifier requried"
		exit 198 
	} 
	if "`mlo'"=="" { 
		di in red "from specifier required"
		exit 198 
	}
	tempname flo fhi val lo hi hold ftry ftry2 tol ltol
	if "`mtol'"=="" & "`mltol'"=="" {
		scalar `tol' = 1e-6
		scalar `ltol' = .
	}
	else if "`mltol'"=="" {
		scalar `tol' = `mtol'
		scalar `ltol'= .
	}
	else if "`mtol'"=="" {
		scalar `tol' = .
		scalar `ltol' =`mltol'
	}
	else {
		scalar `tol' = `mtol'
		scalar `ltol' = `mltol'
	}
	scalar `val' = `mval' 
	scalar `lo' = `mlo' 
	scalar `hi' = `mhi' 
	scalar `ltol' = max(`ltol'*abs(`val'),`ltol')
        local maxiter 64

	if `lo'>`hi' { 
		scalar `hold' = `lo'
		scalar `lo' = `hi'
		scalar `hi' = `hold'
	}

	local ttl : di "Find `progttl'=f() == c=" `val' ", |X_error|<" `tol' /*
		*/ ", |f()-c|<" `ltol'
	local skip = int((72-length("`ttl'"))/2)
	di in gr _n _skip(`skip') "`ttl'" _n

	`callout' quietly `exp1' `lo' `exp2'
	`callout' scalar `flo' = ${`type'} `rtnval' - `val'
	`inline'  scalar `flo' = (`exp1' `lo' `exp2') - `val'

	`callout' quietly `exp1' `hi' `exp2'
	`callout' scalar `fhi' = ${`type'} `rtnval' - `val'
	`inline'  scalar `fhi' = (`exp1' `hi' `exp2') - `val' 

	di in gr /*
		*/ _col(17) "lower" _col(39) "upper" _col(59) "weighted" _n /*
		*/ "iteration" _col(17) "bound" _col(27) "f()-c" /*
		*/ _col(39) "bound" _col(49) "f()-c" /*
		*/ _col(59) "midpoint" _col(72) "f()-c" _n /*
		*/ _dup(76) "-"


	/* verify we bound the solution */
	if `flo'*`fhi'>0 | `flo'*`fhi'==. { 
		di in gr %8.0f 0 ".   " %9.0g `lo' " " %9.0g `flo' /*
			*/ "   " %9.0g `hi' " " %9.0g `fhi'
		display in red "range does not bound solution"
		exit 409
	}
	if `flo'==0 {
		scalar `fhi' = `flo'
		scalar `hi' = `lo'
	}
	else if `fhi'==0 {
		scalar `flo' = `fhi'
		scalar `lo' = `hi'
	}

	local try = (`lo'+`hi')/2
	`callout' quietly `exp1' `try' `exp2'
	`callout' scalar `ftry' = ${`type'} `rtnval' - `val'
	`inline'  scalar `ftry' = (`exp1' `try' `exp2') - `val' 
	local try = `try'+(`try'-`lo')*sign(`flo'-`fhi')*`ftry' / /*
			*/	sqrt(`ftry'^2 - `flo'*`fhi')
	`callout' quietly `exp1' `try' `exp2'
	`callout' scalar `ftry' = ${`type'} `rtnval' - `val'
	`inline'  scalar `ftry' = (`exp1' `try' `exp2') - `val' 
	local iter 0 
	while `ftry'~=0 &( (abs(`hi'-`lo')>2*`tol' | abs(`ftry')>`ltol')) { 
		di in gr %8.0f `iter' ".   " %9.0g `lo' " " %9.0g `flo' /*
			*/ "   " %9.0g `hi' " " %9.0g `fhi' "    " /*
			*/ in ye %9.0g `try' " " %9.0g `ftry'
		if `flo'*`ftry'<0 { 
			scalar `fhi' = `ftry'
			scalar `hi'=`try'
		}
		else {
			scalar `flo' = `ftry'
			scalar `lo' = `try'
		}
		local try = (`lo'+`hi')/2
		`callout' quietly `exp1' `try' `exp2'
		`callout' scalar `ftry' = ${`type'} `rtnval' - `val'
		`inline'  scalar `ftry' = (`exp1' `try' `exp2') - `val'
		local try = `try'+(`try'-`lo')*sign(`flo'-`fhi')*`ftry' / /*
			*/	sqrt(`ftry'^2 - `flo'*`fhi')
		`callout' quietly `exp1' `try' `exp2'
		`callout' scalar `ftry' = ${`type'} `rtnval' - `val'
		`inline'  scalar `ftry' = (`exp1' `try' `exp2') - `val' 
		local iter = `iter' + 1
		if `iter' > `maxiter' {
			global S_1 = `try'
			di in red "Convergence not achieved"
			exit 430
		}
	}
	di in gr %8.0f `iter' ".   " %9.0g `lo' " " %9.0g `flo' /*
		*/ "   " %9.0g `hi' " " %9.0g `fhi' "    " /*
		*/ in ye %9.0g `try' " " %9.0g `ftry'
	di in ye _col(57) "|X_error|<=" %9.0g max(abs(`hi'-`try'), /*
		*/ abs(`try'-`lo'))
	global S_1 = `try'
end
exit

ridder chiprob(2,X)=.05 from 3 to 8 tol 1e-6

-------------------
Note:
        specify:
                rtnval:         S_1             to mean $S_1
                rtnval:         "exp <exp>"     to mean <exp>

        Logic:
                $exp <exp> will evaluate to <exp>


    Find  chiprob(2,X)=f() == c=.05, |X_error|<1.000e-06, |f()-c|<.

         1         2         3         4         5         6         7         8
12345678901234567890123456789012345678901234567890123456789012345678901234567890
                lower                 upper               weighted
iteration       bound     f()-c       bound     f()-c     midpoint     f()-c
----------------------------------------------------------------------------
123456789   123456789 123456789   123456789 123456789    123456789 123456789
                                                        |X_error|<= 1.78e-15
